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ABSTRACT 

We present a powerful and easy-to-implement algorithm for solving constrained op¬ 
timization problems that involve Li/total-variation regularization terms, and both 
equality and inequality constraints. We discuss the relationship of our method to 
earlier works of Goldstein and Osher (2009) and Chartrand and Wohlberg (2010), and 


demonstrate that our approach is a combination of the augmented Lagrangian method 
with splitting and model projection. We test the method on a geomechanical problem 
and invert highly compartmentalized pressure change from noisy surface uplift obser¬ 
vations. We conclude the paper with a discussion of possible extension to a wide class 
of regularized optimization problems with bound and equality constraints. 


INTRODUCTION 

The primary focus of this work is a class of least-squares fitting problems with a total- 
variation (TV) regularization and bound model constraints: 

ll|Vm|||i + ^||F(m)-d||^ ^ min, 

^ (ij 

nil < ni < m 2 . 


In (Q we seek a model vector m such that forward-modeled data F(m) match observed 
data d in the least squares sense, while imposing blockiness-promoting total-variation (TV) 
regularization (Rudin et al. 1992) and lower (mi) and upper (m 2 ) model bounds. Rather 


than using a regularization parameter as a coefficient of the regularization term, we use a 
data-fitting weight a. TV regularization (also know as the Rudin-Osher-Fatemi, or ROF, 
model) acts as a form of model styling that helps to preserve sharp contrasts and boundaries 
in the model, even when spectral content of input data has limited resolution. Examples 
of successful geophysical application of unconstrained TV-regularized optimization are in¬ 
cluded in (Maharramov and Biondi 2015; Maharramov et al., 2015 Ma et al. 2015a|b). The 


regularization provided by bounded total-variation sometimes produces sufficient smoothing 
side-effect on the inverted model that obviates explicit bound constraints. However, many 
applications still require the imposition of additional constraints regardless of the regular¬ 
ization. For example, reservoir pore-pressure inversion problems often come with a priori 
bounds on the estimated pore pressure change, such as the pore pressure change being non¬ 
negative as a result of fluid injection (lower bound) or never exceeding a hydraulic fracturing 
pressure (upper bound). An example of such inversion for an unconventional reservoir from 


field tilt measurements is provided by Maharramov and Zoback (2015). TV regularization 


is a key tool in imaging and de-noising applications (Rudin et al. 1992 Chambolle and 
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Lions, 1997 Goldstein and Osher| 2009; Chartrand and Wohlberg, 2010) and require an 


efficient mechanism for including a priori model constraints that can significantly reduce 
model space ( Chartrand and Wohlberg|2010 ). While barrier or penalty function methods, 
such as nonlinear interior-point methods ( [Nocedal and Wright , 2006), can be used to tackle 
the general constrained formulation (|^, the presence of a non-differentiable Li-norm total- 
variation term and non-quadratic penalty terms pose considerable challenges to practical 
implementation. A log-barrier function such as 


const X E log 


1712 — m} 


+ log 


m* — m\ 


( 2 ) 


2=1 


where n is the model space dimension, can be added to the right-hand side of the objective 
function to keep solution iterates away from the rectangular bounds. However, this adds 
a non-quadratic term to the objective function. For large-scale inversion problems with 
n > 10^ (such as typical in geophysical applications) often only iterative gradient-based 


solution techniques like the nonlinear conjugate gradients (Nocedal and Wright, 2006) are 


available, and adding non-quadratic terms may significantly affect convergence properties. 
Note that this is in addition to the challenges associated with handling the non-differentiable 
TV-regularization term. 


Chartrand and Wohlberg (2010) used a splitting approach to decouple the TV-regularized 


problem from enforcing the constraints. In our approach, we perform three-way splitting of 
problem ([T|) into a smooth optimization, gradient thresholding and projection steps using 
the Alternating Direction Method of Multipliers (ADMM) (Boyd et al., 2010j). For uncon¬ 


strained TV-regularized problems this approach is equivalent to the split-Bregman method 


of Goldstein and Osher (2009). However, we integrate the projection step associated with 


enforcing the bound constraints into the TV-minimization loop and avoid unnecessary it¬ 


erations in the minimization of a proximal term (Parikh and Boyd, 2013) associated with 
the projection. 


METHOD 


First, we recast the TV-regularization part of ( I|) as a constrained optimization problem 
following the approach of Goldstein and Osher (2009). We introduce an auxiliary variable 
X and operator $ : m —)• x such that for isotropic TV regularization we have a vector of 
the model-space dimension 


$(m) = y (Vj:m)^ -F (Vym)^, 

and for anisotropic regularization a vector twice the model-space dimension 

Vajin 


Vyin 


$(m) = 

Problem ([^ can now be reformulated with an additional equality constraint; 


( 3 ) 


( 4 ) 


a , 


l|x||i + 2 ^ 

X = $(m), 
mi < m < m 2 . 


mm, 


( 5 ) 
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Problem ([^ is still a bound-constrained problem. Introducing the projection operator 

n(m) = max{min{m, m 2 }, mi}, (6) 

where min and max are applied component-wise, we reduce to a fully equality-constrained 
formulation: 


l|x||i + -||F(m)-d ||2 
X = $(m), 
m = y. 


mm, 


( 7 ) 


y = n(m). 


Following the augmented Lagrangian recipe for 0 while assuming the last equality con¬ 
straint still enforced explicitly, we obtain a sequence of problems (Nocedal and Wrightf 
200^ 


X 


fc + l _ 


a , 


, m 


= argmin ||x||i -|- — ||F(m) — d|| 2 -|- 


A (5 

-||x - $(m)||| -/^fc^(x - $(m))-h-||m-y|||-i/fc^(m-y) 


mm, 


Mfc+i = - A - $(m''+^) 




( 8 ) 




m^+1 - y 


A: = 0,1,2,... 


Coefficients A and 5 are any positive constants above certain problem-specific “threshold” 
values (Nocedal and Wright [2006[ ), and can be selected experimentally. Vectors /x^ and i^k 
are vectors of multipliers that converge to the set of Lagrange multipliers for the first two 
equality constraints of problem ([^. At each step, ([^ solves an Li-regularized problem with 
respect to the combined model vector (x, m) . Introducing new scaled multiplier vectors 


= 


f, c‘ = J,k = 0,l,2, 


( 9 ) 


a little algebra shows that 
(x^+^ 


is equivalent to 


,m 


k+l\ _ 


II II ^11 / \ 1 I I 2 

= argmin ||x||i -|- —||F(m) — d|| 2 -|- 


^||x - $(m)-b^||i + ^||m - y-c^lli 
bfc+i ^ b^ + $(m^+i)-x^+\ 

^ + /c = 0,1,2, 


mm, 


( 10 ) 


Here we used the fact that adding a constant term A/2||b^||2-|-<j/2||c^||2 to the objective 
function obviously does not change the minimizing solution. 

Problem ([^ can be solved by iteratively projecting the current model vector m onto y, 
then conducting the iterations (10) to convergence, then repeating the process. However, 
presence of the proximal term (l/2||m — yH^ in 


due to the constraint m = y means that 
a very accurate solution of at early iterations is wasteful and unnecessary. We instead 
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carry out a single iteration of ( 10 ) followed by the model projection: 

|2 


(X^ 


k-\-l _ 


a , 


= argmin ||x||i + —||F(m) — d|| 2 + 


^||x - $(m) - b^lll +-||m - y"-c "||2 


-.k\\2 


mm, 


k-\-l _ 


= b'' + $(m 


k+i)_^k+i 


( 11 ) 


c^+i = c^ + y^-m^+\ 


,k+l 


y""^ = n( ) = maxlminlm^"^ ,1112}, mi}, /c = 0 , 1 , 2 ,... 

The iterative process 01 still requires soling an Li-regularized problem. However, the 
Li-norm term now involves only the vector x. Therefore, we apply splitting, minimizing 


. ,,|2 A, 


|x||i + ^||F(m)-d||^ +^||x - $(m)-b'=||i + -||m - y"^ - 


( 12 ) 


alternately with respect to m and x in an inner loop of A^i > 1 cycles. Because the proximal 
constraint m = y renders good fitting accuracy at early stages unnecessary, A^i can be small. 


Further we note that the minimization of (12) with respect to x (in a splitting step with m 


fixed) is given trivially by the “shrinkage” operator (Goldstein and Osher, 2009): 


x^^^ = shrink $(m) + b*^, — } , 


where 


shrink{x, 7 } = -j—r max (|x| — 7 , 0 ) 



(13) 

7,0), 

(14) 


and is effectively thresholding the model gradient. Our algorithm can be described by the 
following 5 steps: 


1 Initialization 

= starting guess, 

x° = 0, 

y^ = max{min{m*^, m 2 }, mil, 
b° = 0, 
c° = 0, 

2 Outer loop. Repeat steps 3-5 for A; = 0,1, 2,... 


(15) 


3 Inner loop. Iterate (16) A^i > 1 times 

A 


m'"'" = argmin — $(m) — b *^||2 -|- ^||F(m) — d|||-|- 

^11 h r ii 9 

-||m-y'' 


(16) 


X 


/c+1 


= shrink < $(m' 




+ b^7 


x^ = x^+\ 
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4 Update the multipliers and project the model onto the bounding rectangle: 

bfc+i ^ 

^fc+i ^ c'=+y^-m^+\ (17) 

= max{min{m^"''^, m 2 }, mi}. 


5 Terminate if the target accuracy is reached 

< target accuracy. 


| m *^+^ — m ^||2 


m'' 


or go back to step 2 otherwise. 


(18) 


Optimizing (16) with respect to m is in itself a large-scale optimization problem, nonlinear 
for a nonlinear modeling operator F. Solving the optimization problem (16) exactly is 
unnecessary because for small k (i.e., at early stages of the inversion) vector y^ is not the 
true model, vector x^ is far from the true model gradient, and the multipliers b*^,x^ could 
be far from scaled Lagrange multipliers, s. Therefore, for large-scale problems only a few 
steps of an iterative method like conjugate gradients need be carried out. As the solution 
converges to the true solution and critical sharp contrasts in the model are identified, an 
iterative solver can begin to take advantage of the objective function curvature information 
collected at previous iterations of the outer loop, potentially leading to a significantly faster 
convergence. Optimal strategies for spanning iterations of nonlinear conjugate gradients 
across iterations of the outer loop of our algorithm are the subject of an upcoming report. 


RESULTS 


We demonstrate our method with a test problem that simulates vertical surface uplift in re¬ 
sponse to distributed dilatational sources, mathematically equivalent to surface deformation 


due to pore pressure change (Segall 2010). Our modeling operator is defined as 


F(m) = u{x), u{x) = / 

Jo 


D^m{^)d^ 


0 (7l2 + (a;_^)2)3/2’ 


(19) 


where we assume that m = m(^),^ G [0;^] is a relative pore pressure change along a 
linear segment [0, A] of a reservoir at a constant depth D, and u = u{x),x G [0, A], within 


a proportionality factor determined by poroelastic medium properties (Maharramov and 


Zoback 2015), is the induced vertical uplift on the surface. For demonstration purposes we 


consider a one-dimensional model but the results trivially extend to realistic reservoir and 


surface geometries. Operator (19) is a smoothing operator, and recovering sharp pressure 
contrasts e.g. due to permeability barriers requires model “styling” or regularization such 
as blockiness-promoting ROF model. As a true model we used a highly compartmentalized 
pressure model of Figure |lb[ In our experiments, we set D = 100m A = 2km, and discretized 
both the model and data space using a 200-point uniform grid. Random Gaussian noise with 
a = 15% of the maximum clean data amplitude was added to the clean forward-modeled 


data to produce the noisy observations shown in Figure la 


The result of a TV-regularized unconstrained inversion is shown in Figure 2a against 
the true model and a Tikhonov-regularized inversion. This result was obtained using the 
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True and noisy uplift observations 



X (km) 


oc 


True model 


X (km) 


a b 

Figure 1: (a) True and noisy uplift observations. Random Gaussian noise with a = 15% of 
maximum clean data amplitude was added to the clean data, (b) True model exhibibits a 
highly compartmentalized “blocky” behavior. 


above algorithm by setting <5 = 0 (no bound constraints) and using the values of a = 1 
and A = 2. The TV-regularized result captures the compartmentalized picture of pressure 
distribution better than the highly smoothed Tikhonov regularization result. However, 
due to absence of bound constraints, lower pressure bounds are not honored, resulting in 
negative pressure areas that are not present in the true model. The result of running our 
bound-constrained TV-regularization algorithm is shown in Figure 2b The imposition of 


bound constraints not only removed the negative relative pressure areas, but also removed 
the pressure spike at about x ~ 1km in the unconstrained inversion of Figure that 
apparently had resulted from compensating negative pressures. In both the constrained 
and the unconstrained runs we conducted 1000 outer loop iterations with 2 inner loops 
cycles. However, the algorithm converged quickly, with only a few initial iterates outside a 
tight neighborhood of the final curve, as shown in Figure 3b Finally, we note that many 


practical implementations of bound constraints often resort to a simplistic way of enforcing 
the constraints: the inverted model is projected onto the bounding rectangle either once 
after applying a direct unconstrained solver, or at each iteration of an unconstrained solver. 
In this case variable y and the associated quadratic regularization term are not introduced 
into the objective function. This may result in a violation of the KKT optimality conditions 
where the bound constraints are active (Nocedal and Wright, 2006| 
by the blue plot in Figure While the bound constraints are honored 


and is demonstrated 
the solution is 


both qualitatively and quantitatively far from optimal. 


CONCLUSIONS AND PERSPECTIVES 

Our algorithm combines the advantages of the blockiness-promoting and edge-preserving 
ROF model with the ability to impose bound constraints. The splitting mechanism used 
for enforcing the bound constraints is naturally integrated into the split-Bregman iterations 
and results in no extra computational cost. The method was able to resolve compartmen- 
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talized subsurface pressure changes from noisy surface uplift observations despite the highly 
diffusive nature of the underlying deformation process. The method can be implemented 


Inversion results Inversion results 




a b 

Figure 2: (a) Unconstrained TV-regularized inversion. The algorithm tries to fit the data by 
allowing negative relative pressure changes, (b) Bound constrained TV-regularized result. 
Note that enforcing lower bounds resulted in a more accurate shape matching of the true 
model. 

around any large-scale nonlinear solver such as conjugate gradients or quasi-Newton meth¬ 
ods. Additional equality and inequality constraints can be incorporated into the algorithm 
using the general ADMM framework. 
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Inversion results Convergence of TV-regulazed models 




a b 

Figure 3: (a) Direct imposition of the bound constraints at each iteration of the uncon¬ 
strained solver resulted in a qualitatively and quantitatively wrong inversion, (b) Con¬ 
vergence of TV-regularized inverted models with bound constraints. The method quickly 
resolves both sharp contrasts and active bounds as only a few initial curves out of 1000 
iterates lie outside a small neighborhood of the final curve. 
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